package scu.maqiang.numeric;

import java.util.Arrays;
import java.util.stream.DoubleStream;

/**
 * 统计类, 实现基本的统计计算
 *
 * @author 马强
 */
public class Statistics {

//	public static double[] t_0975 = {0.0, 12.7062, 4.3027, 3.1824, 2.7764, 2.5706, 
//			                        2.4469, 2.3646, 2.3060, 2.2622, 2.2281,
//			                        2.2010, 2.1788, 2.1604, 2.1448, 2.1315, 
//			                        2.1199, 2.1098, 2.1009, 2.0930, 2.0860, 
//			                        2.0796, 2.0739, 2.0687, 2.0639, 2.0595, 
//			                        2.0555, 2.0518, 2.0484, 2.0452, 2.0423, 
//			                        2.0395, 2.0369, 2.0345, 2.0322, 2.0301,
//			                        0.0281, 2.0262, 2.0244, 2.0227, 2.0211,
//			                        2.0195, 20.181, 2.0167, 2.0154, 2.0141};

    public static double[] t_0975 = {0.0, 12.706204736174696, 4.302652729749463, 3.182446305283710, 2.776445105197793, 2.570581835636315,
            2.446911851144969, 2.364624251592784, 2.306004135204166, 2.262157162798204, 2.228138851986273,
            2.200985160091638, 2.178812829667226, 2.160368656462786, 2.144786687917803, 2.131449545559775,
            2.119905299221255, 2.109815577833316, 2.100922040241036, 2.093024054408310, 2.085963447265863,
            2.079613844727680, 2.073873067904024, 2.068657610419049, 2.063898561628026, 2.059538552753295,
            2.055529438642873, 2.051830516480287, 2.048407141795244, 2.045229642132705, 2.042272456301238,
            2.039513446396409, 2.036933343460100, 2.034515297449339, 2.032244509317719, 2.030107928250342,
            2.028094000980449, 2.026192463029110, 2.024394163911969, 2.022690920036762, 2.021075390306272,
            2.019540970441375, 2.018081702818444, 2.016692199227825, 2.015367574443762, 2.014103388880851,
            2.012895598919424, 2.011740513729769, 2.010634757624234, 2.009575237129235, 2.008559112100758};

    public static double[] t_0950 = {0.0, 6.313751514675031, 2.919985580353724, 2.353363434801823, 2.131846786326649, 2.015048373333023,
            1.943180280515303, 1.894578605090007, 1.859548037530899, 1.833112932656237, 1.812461122811674,
            1.795884818704040, 1.782287555649314, 1.770933395986864, 1.761310135774879, 1.753050355692573,
            1.745883676276249, 1.739606726075073, 1.734063606617537, 1.729132811521370, 1.724718242920788,
            1.720742902811880, 1.717144374380241, 1.713871527747050, 1.710882079909429, 1.708140761251896,
            1.705617919759272, 1.703288445722127, 1.701130934265929, 1.699127026533497, 1.697260886593960,
            1.695518782545866, 1.693888748383711, 1.692360309030344, 1.690924255186855, 1.689572457780266,
            1.688297714116813, 1.687093619596263, 1.685954460166749, 1.684875121711213, 1.683851013335661,
            1.682878002132711, 1.681952357467542, 1.681070703202515, 1.680229976572132, 1.679427392652322,
            1.678660413556916, 1.677926721641809, 1.677224196124314, 1.676550892616896, 1.675905025163113};

    public static double Phi_0975 = 1.96;
    public static double Phi_0950 = 1.645;

    public static double[] Chi2_0975 = {0.0, 5.023886187314888, 7.377758908227871, 9.348403604496145, 11.143286781877794, 12.832501994030023,
            14.449375335447918, 16.012764274629323, 17.534546139484647, 19.022767798641635, 20.483177350807395,
            21.920049261021205, 23.336664158645341, 24.735604884931568, 26.118948045037367, 27.488392863442972,
            28.845350723404753, 30.191009121639805, 31.526378440386623, 32.852326861729701, 34.169606902838339,
            35.478875905727250, 36.780712084035549, 38.075627250355794, 39.364077026603908, 40.646469120275192,
            41.923170096353914, 43.194510966156038, 44.460791836317760, 45.722285804174540, 46.979242243671166,
            48.231889594451943, 49.480437742971688, 50.725080066281237, 51.965995195121913, 53.203348542056489,
            54.437293631813219, 55.667973264261100, 56.895520535055986, 58.120059734686279, 59.341707143171206,
            60.560571734843755, 61.776755805349204, 62.990355531101976, 64.201461469886752, 65.410159009999575,
            66.616528774250455, 67.820646984252463, 69.022585789666081, 70.222413566434525, 71.420195187506422};

    public static double[] Chi2_0025 = {0.0, 0.000982069117175, 0.050635615968580, 0.215795282623898, 0.484418557087930, 0.831211613486658,
            1.237344245791202, 1.689869180677354, 2.179730747252649, 2.700389499980358, 3.246972780236841,
            3.815748252236099, 4.403788506981701, 5.008750511810332, 5.628726103039733, 6.262137795043251,
            6.907664353497004, 7.564186449577570, 8.230746194756662, 8.906516481987975, 9.590777392264869,
            10.282897782522863, 10.982320734473674, 11.688551922452435, 12.401150217444432, 13.119720024937781,
            13.843904982007606, 14.573382730821706, 15.307860552601193, 16.047071695364885, 16.790772265566623,
            17.538738581475489, 18.290764907283041, 19.046661503175120, 19.806252939214573, 20.569376630744969,
            21.335881560799059, 22.105627161169515, 22.878482328733469, 23.654324557593018, 24.433039170807891,
            25.214518638112501, 25.998661968152369, 26.785374165536332, 27.574565744459218, 28.366152291859844,
            29.160054074089363, 29.956195681912103, 30.754505709372932, 31.554916462667123, 32.357363695658641};

    /**
     * 计算统计量的平均值
     *
     * @param x 统计值
     * @return 平均值
     */
    public static double mean(double... x) {
        return DoubleStream.of(x).average().orElse(Double.NaN);
    }

    /**
     * 计算统计量的样本标准差
     *
     * @param x 统计值
     * @return 样本标准差
     */
    public static double std(double... x) {
        double mu = mean(x);
        int n = x.length;
        return Math.sqrt(DoubleStream.of(x).map(ele -> (ele - mu) * (ele - mu)).sum() / (n - 1));
    }

    public static double[] muCI_t(double mean, double std, int n) {
        double temp = t_0975[n - 1] * std / Math.sqrt(n);
        return new double[]{mean - temp, mean + temp};
    }

    public static double[] muCI_normal(double mean, double std, int n) {
        double temp = Phi_0975 * std / Math.sqrt(n);
        return new double[]{mean - temp, mean + temp};
    }

    public static double[] sigmaCI(double std, int n) {
        return new double[]{std * Math.sqrt((n - 1) / Chi2_0975[n - 1]), std * Math.sqrt((n - 1) / Chi2_0025[n - 1])};
    }

    public static double[] sigmaCI(double[] x) {
        double myStd = std(x);
        System.out.println(myStd);
        return sigmaCI(myStd, x.length);
    }


    public static double[] muCI(double[] x, Double STD) {
        double mean = mean(x);
        System.out.println(mean);
        if (STD == null) {
            double sampleStd = std(x);
            return muCI_t(mean, sampleStd, x.length);
        } else {
            return muCI_normal(mean, STD, x.length);
        }
    }

    public static double GammaFunc(double x) {
        int i;
        double y, t, s, u;
        double[] a = {0.0000677106, -0.0003442342, 0.0015397681, -0.0024467480, 0.0109736958, -0.0002109075, 0.0742379071, 0.0815782188, 0.4118402518, 0.4227843370, 1.0};
        if (x <= 0.0) {
            System.out.println("Error * * x <= 0!");
            return (-1.0);
        }
        y = x;
        if (y < 1.0) {
            t = 1.0 / (y * (y + 1.0));
            y += 2.0;
        } else if (y <= 2.0) {
            t = 1.0 / y;
            y += 1.0;
        } else if (y <= 3.0) {
            t = 1.0;
        } else {
            t = 1.0;
            while (y > 3.0) {
                y -= 1.0;
                t *= y;
            }
        }
        s = a[0];
        u = y - 2.0;
        for (i = 1; i <= 10; i++) {
            s = s * u + a[i];
        }
        s *= t;
        return (s);
    }

    public static double BetaFunc(double x, double y) {
        if (x <= 0.0 || y <= 0.0) {
            System.out.println("Error * * x <= 0 or y <= 0.0!");
            return (-1.0);
        }
        return GammaFunc(x) * GammaFunc(y) / GammaFunc(x + y);
    }

    public static double Bessel_1(int n, double x) {
        int i, m;
        double t, y, z, p, q, s, b0, b1;
        p = 0.0;
        double[] a = {57568490574.0, -13362590354.0, 651619640.7, -11214424.18, 77392.33017, -184.9052456};
        double[] b = {57568490411.0, 1029532985.0, 9494680.718, 59272.64853, 267.8532712, 1.0};
        double[] c = {72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.4826, -30.16036606};
        double[] d = {144725228443.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0};
        double[] e = {1.0, -0.1098628627e-02, 0.2734510407e-04, -0.2073370639e-05, 0.2093887211e-06};
        double[] f = {-0.1562499995e-01, 0.1430488765e-03, -0.6911147651e-05, 0.7621095161e-06, -0.934935152e-07};
        double[] g = {1.0, 0.183105e-02, -0.3516396496e-04, 0.2457520174e-05, -0.240337019e-06};
        double[] h = {0.4687499995e-01, -0.2002690873e-03, 0.8449199096e-05, -0.88228987e-06, 0.105787412e-06};
        t = Math.abs(x);
        if (n < 0) n = -n;
        if (n != 1) {
            if (t < 8.0) {
                y = t * t;
                p = a[5];
                q = b[5];
                for (i = 4; i >= 0; i--) {
                    p = p * y + a[i];
                    q = q * y + b[i];
                }
                p = p / q;
            } else {
                z = 8.0 / t;
                y = z * z;
                p = e[4];
                q = f[4];
                for (i = 3; i >= 0; i--) {
                    p = p * y + e[i];
                    q = q * y + f[i];
                }
                s = t - 0.785398164;
                p = p * Math.cos(s) - z * q * Math.sin(s);
                p = p * Math.sqrt(0.636619772 / t);
            }
        }
        if (n == 0) return (p);
        b0 = p;
        if (t < 8.0) {
            y = t * t;
            p = c[5];
            q = d[5];
            for (i = 4; i >= 0; i--) {
                p = p * y + c[i];
                q = q * y + d[i];
            }
            p = x * p / q;
        } else {
            z = 8.0 / t;
            y = z * z;
            p = g[4];
            q = h[4];
            for (i = 3; i >= 0; i--) {
                p = p * y + g[i];
                q = q * y + h[i];
            }
            s = t - 2.356194491;
            p = p * Math.cos(s) - z * q * Math.sin(s);
            p = p * x * Math.sqrt(0.636619772 / t) / t;
        }
        if (n == 1) return (p);
        b1 = p;
        if (x == 0.0) return (0.0);
        s = 2.0 / t;
        if (t > 1.0 * n) {
            if (x < 0.0) b1 = -b1;
            for (i = 1; i <= n - 1; i++) {
                p = s * i * b1 - b0;
                b0 = b1;
                b1 = p;
            }
        } else {
            m = (n + (int) Math.sqrt(40.0 * n)) / 2;
            m = 2 * m;
            p = 0.0;
            q = 0.0;
            b0 = 1.0;
            b1 = 0.0;
            for (i = m - 1; i >= 0; i--) {
                t = s * (i + 1) * b0 - b1;
                b1 = b0;
                b0 = t;
                if (Math.abs(b0) > 1.0e+10) {
                    b0 = b0 * 1.0e-10;
                    b1 = b1 * 1.0e-10;
                    p = p * 1.0e-10;
                    q = q * 1.0e-10;
                }
                if ((i + 2) % 2 == 0) q = q + b0;
                if ((i + 1) == n) p = b1;
            }
            q = 2.0 * q - b0;
            p = p / q;
        }
        if ((x < 0.0) && (n % 2 == 1)) p = -p;
        return (p);
    }

    public static void main(String[] args) {
        double[] y = {0.39, 0.40, 0.50, 0.39, 0.4, 0.36, 0.37, 0.42, 0.4, 0.37};
        System.out.println("Mean: " + mean(y));
        System.out.println("Std: " + std(y));
        y = new double[]{14.5, 17.4, 16.7, 15.2, 17.3, 14.7, 17.1, 16.0, 14.8, 13.7};
        System.out.println("Mean: " + mean(y));
        System.out.println("Std: " + std(y));

        double[] CI = muCI_t(0.95, 0.4, 25);
        System.out.println("Confidential Interval: " + Arrays.toString(CI));

        double[] sample = {0.19, 0.24, 0.04, 0.08, 0.20, 0.12, 0.31, 0.29, 0.13, 0.07};
        CI = muCI(sample, 0.16);


        System.out.println("Confidential Interval: " + Arrays.toString(CI));
//		double[] sigmaCI = sigmaCI(sample);
//		System.out.println("Confidential Interval: " + Arrays.toString(sigmaCI));
//		
        double[] sample2
                = {578, 572, 570, 568, 572, 570, 570, 596, 584, 572};
        CI = muCI(sample2, 5.0);
        System.out.println("Confidential Interval: " + Arrays.toString(CI));
        CI = muCI(sample2, null);
        System.out.println("Confidential Interval: " + Arrays.toString(CI));

        for (int i = 0; i <= 10; i++) {
            double x = 0.5 * i;
            double gx = GammaFunc(x);
            System.out.println("x = " + x + "    gamma(x) = " + gx);
        }
        System.out.println("Beta(1.5, 2.5) = " + BetaFunc(1.5, 2.5));

        for (int n = 0; n <= 5; n++) {
            double x = 0.05;
            for (int i = 1; i <= 4; i++) {
                double bnx = Bessel_1(n, x);
                System.out.println("n=" + n + "   x= " + x + "    J(n, x)= " + bnx);
                x = x * 10.0;
            }
        }
//		
//		double[] sample3 = {10.1, 10.3, 10.4, 10.5, 10.2, 9.7, 9.8, 9.9, 10.0, 10.1, 9.8, 10.2};
//		CI = muCI(sample3, null);
//		System.out.println("Confidential Interval: " + Arrays.toString(CI));
//		sigmaCI = sigmaCI(sample3);
//		System.out.println("Confidential Interval: " + Arrays.toString(CI));

//		double[] sample4 = {11.5, 11.21, 11.05, 11.08, 11.07, 11.1, 11.06, 11.04, 11.11, 11.03};
//		double[] sigmaCI = sigmaCI(sample4);
//		System.out.println("Confidential Interval: " + Arrays.toString(sigmaCI));
    }
}
